#include "ReadDelMapping.h"

const unsigned int BUFFERED_READS_SIZE = 1000000;

CDelMappingsQ::CDelMappingsQ(void)
{
    this->clear();
}

CDelMappingsQ::~CDelMappingsQ(void)
{
}

void CDelMappingsQ::clear(void)
{
    this->minDiff = INT_MAX;
    this->minDisScore = INT_MAX;
    this->delAlignments.clear();
    this->bAmbiguousDelMapping = false;
}

void CDelMappingsQ::push_back(CMappingResult &m)
{
    if ((int)this->delAlignments.size() < this->MAX_ALIGNMENT_NO) {
        this->delAlignments.push_back(m);
    } else { // replace a worse record with current better record
        vector<CMappingResult>::iterator it = this->delAlignments.begin();
        for (; it != this->delAlignments.end(); it++) {
            if (m.mismatchScore < it->mismatchScore) {
                //*it = m;
                // upgrade is update the new alignment of a same read.
                it->upgrade(m);
                break;
            }
        }
    }
}

void CDelMappingsQ::update(CMappingResult &m, bool bUniqueDelMapOnly, bool bAllDelMap)
{
    if (bAllDelMap) {
        this->push_back(m);
        if (this->minDisScore > m.mismatchScore) {
            this->minDisScore = m.mismatchScore;
            this->minDiff = m.uiDiff;
            this->bAmbiguousDelMapping = false;
        } else if (this->minDisScore == m.mismatchScore) {
            this->bAmbiguousDelMapping = true;
        }
    } else {
        if (this->minDisScore > m.mismatchScore) {
            this->clear();
            this->minDisScore = m.mismatchScore;
            this->minDiff = m.uiDiff;
            this->bAmbiguousDelMapping = false;
            this->delAlignments.push_back(m);
        } else if (this->minDisScore == m.mismatchScore) {
            this->bAmbiguousDelMapping = true;
            if (!bUniqueDelMapOnly) {
                this->push_back(m);
            }
        }
    }
}

CReadsDelMapping::CReadsDelMapping(void)
{
    this->initialization();
}

CReadsDelMapping::~CReadsDelMapping(void)
{
}

CReadsDelMapping::CReadsDelMapping(const MappingOpts P) : CReadsMapping( P )
{
    this->initialization();
}

void CReadsDelMapping::initialization(void)
{
    this->noDelMapping = 0;
    this->noOfAmbiguousDelMapping = 0;
    this->noOfDelSepLess = 0;
    this->noOfDelSepMore = 0;
    this->noOfDelSepMoreAndLess = 0;
}

int CReadsDelMapping::locateDelAlignments
(const CGenome_Index_TableQ& table,CLongReadInBits & read,  const char* readName, const char* pQScore)
{
    // 1. Search head and tail
    alignmentsQ[0].read = read.head();
    alignmentsQ[1].read = read.tail();

    table.queryReadBases(alignmentsQ[0].read, alignmentsQ[0], true, true);
    table.queryReadBases(alignmentsQ[0].read, alignmentsQ[0], false, false);
    table.queryReadBases(alignmentsQ[1].read, alignmentsQ[1], true, true);
    table.queryReadBases(alignmentsQ[1].read, alignmentsQ[1], false, false);

    // 2. Get reads, tag and quality score
    CMappingResult m(alignmentsQ[0], read.uiSeqLength);
    myStrCpy(m.QNAME, readName, FILENAME_MAX);
    read.decode(m.caRead, READ_BUFFER_SIZE);
    if (pQScore != NULL) {
        if (read.uiSeqLength < MAX_READ_LENGTH) {
            trQScores(read.uiSeqLength, 0, pQScore, m.rawScores);
            trQScores(read.uiSeqLength, SolexaScoreEncodingShift, m.rawScores, m.QScores);
        }
    }
    this->examDelAlignments4aRead(table, m);

    return(0);
}

// TODO pass in special data structure so it can checked
int CReadsDelMapping::getDelAlignments(const CGenome_Index_TableQ& table, CLongReadsSet& readSet)
{
    if (table.bMapReadInColors) {
        cout << "Currently Clip doesn't support mapping SOLiD reads" << endl;
        return(-1);
    }
    // Flag that set the alignment is ambiguous or not
    getReadsFileFormat(readSet.InputFile, opt.readsFileFormat);
    unsigned int uiReadLength = readSet.uiRead_Length;
    string seedStr = seedSymbol(table.chosenSeedId);
    printf("Mapping %s (%d-bp reads) with %s seed.\n", \
           readSet.InputFile, uiReadLength, seedStr.c_str());

    this->initializeStatsCounter();
    this->setUpIO4Aligment( readSet.InputFile, table);
    bool hasQScore = (readSet.cFileType == 'q' && !readSet.bIgnoreQScores);
    if (hasQScore) {
        int AVG_SCORE_PER_BASE = 20;
        this->opt.subDiffThreshold *= AVG_SCORE_PER_BASE;
    }
    while (readSet.get_next_capacity_reads(BUFFERED_READS_SIZE) > 0) {
        int readSetSize = (int)readSet.longReadsInBits.size();
        for (int i = 0; i < readSetSize; i++) {
            const char* pQScore = hasQScore ? readSet.pQualScores->qScores(i) : NULL;
            // Query read
            CLongReadInBits& read = readSet.longReadsInBits.at(i);
            if (i % this->CHECK_POINTS == 0) {
                printf("Mapping no %u reads.\r", iReadCounter + i);
            }
            locateDelAlignments(table, read, readSet.readsId.at(i).id, pQScore);
        }
        this->iBadReadCounter += (int)readSet.uiNo_of_Bad_Reads;
        this->iReadCounter += (int)readSet.readsId.size();
        printf("Mapping no %u reads.\r", iReadCounter);
    }
    this->tearDownIO4Aligment();
    this->iMapCount = this->noDelMapping;
    this->printMappingStats(cout, readSet.InputFile, opt.subDiffThreshold);
    //this->printInDelAlignStats(cout, string(readSet.InputFile), opt.subDiffThreshold);
    return(0);
}

// start is the start of ths substring, length is the length of the substring for comparision
int getMismatchScore(const char* ref, const char* seq, const char* scores, int start, int length)
{
    // Note it is still possible to overflow in the ref, seq and scores array
    int score = 0;
    if (ref && seq && (start >= 0) && (length >= 0)) {
        bool bQScoreAvailable = !(scores == NULL || scores[0] == '\0');
        for (int i = start; i < start + length; i ++ ) {
            if (ref[i] == '\0' || seq [i] == '\0') {
                ERR
                break;
            } else if ( ref[i] != seq[i]) {
                if (bQScoreAvailable) {
                    score += (int)(scores[i]);
                } else {
                    score ++;
                }
            }
        }
    } else {
        ERR
    }
    return(score);
}

inline int log4(int num)
{
    int returnV = 0;
    while (num > 0) {
        num /= 4;
        returnV++;
    }
    return(returnV);
}

inline int myChoose(int n, int m)
{
    switch (m) {
    case 0:
        return(1);
    case 1:
        return(n);
    case 2:
        return( n * (n - 1)/ 2 );
    case 3:
        return( n * (n - 1) * (n - 2) / 6 );
    case 4:
        return( n * (n - 1) * (n - 2) * (n - 3) / 24 );
    case 5:
        return( n * (n - 1) * (n - 2) * (n - 3) * (n - 4) / 120 );
    case 6:
        return( n * (n - 1) * (n - 2) * (n - 3) * (n - 4) * (n - 6) / 720 );
    default:
        return(0);
    }
}

inline int getMatchScore(int length, int mismatches)
{
    const int minScoreThreshold = 5;
    const int maxScoreThreshold = 30;
    int matchScore = length - mismatches;
    if (length <= minScoreThreshold ) {
        return(0);
    } else if (matchScore > maxScoreThreshold) {
        matchScore = maxScoreThreshold;
    } else if (mismatches == 0) {
        matchScore = length;
    } else {
        int chooseV = myChoose(length, mismatches);
        if (chooseV) {
            matchScore -= log4(chooseV);
        } else {
            matchScore -= (2 * mismatches);
        }
    }
    if ( matchScore > maxScoreThreshold) {
        matchScore = maxScoreThreshold;
    } else if (matchScore <= 0) {
        return(0);
    }
    return( matchScore);
}

/*
 *  Given a ref, a read, and the possible del alignment specified by the two mapped end in m,
 *  find the best splicing (Del) position of the alignment with mismatch score lower than the threshold misT.
 */
bool CReadsDelMapping::examDelAlignment(CSeqInBits& ref, CMappingResult& m, unsigned int misT)
{
    // Here I implement a non-bitwise operation version which can be improved for speed later
    char refBuf1[MAX_LINE];
    char refBuf2[MAX_LINE];
    if (m.uiReadLength < this->opt.anchorLength) {
        ERR
    }
    unsigned int refL = m.uiRefLength;
    unsigned int refIndexShift;
    if (m.uiGlobalMappedPos >= m.uiPOS) {
        refIndexShift = m.uiGlobalMappedPos - m.uiPOS; // global start index of the ref
    } else {
        LOG_INFO("Info %d: Inconsistent mapped index.\r", WARNING_LOG);
        refIndexShift = max(m.uiGlobalMappedPos, m.uiPOS);
    }
    unsigned int examLength = m.uiReadLength - this->opt.anchorLength;
    unsigned int anchorL = this->opt.anchorLength;
    unsigned int headSubStrStart = m.uiGlobalMappedPos;
    unsigned int headSubStrEnd = m.uiGlobalMappedPos + examLength;
    unsigned int tailSubStrStart = m.uiMPOS + refIndexShift + this->opt.anchorLength - examLength;
    unsigned int tailSubStrEnd = m.uiMPOS + refIndexShift + this->opt.anchorLength;

    if ( headSubStrEnd > refIndexShift + refL ||\
            tailSubStrEnd > refIndexShift + refL ||\
            headSubStrStart < refIndexShift ||\
            headSubStrStart > headSubStrEnd || tailSubStrStart > tailSubStrEnd ||\
            headSubStrStart >= tailSubStrStart) {
        return(false);
    }

    bool bGoodDelMap = false;
    char* headRef = ref.decodeSubStr(headSubStrStart, headSubStrEnd, refBuf1, MAX_LINE);
    char* tailRef = ref.decodeSubStr(tailSubStrStart, tailSubStrEnd, refBuf2, MAX_LINE);
    if (m.strand == '-') {
        // reverse the reference
        reverseComplementKmer(headRef);
        reverseComplementKmer(tailRef);
        swap(headRef, tailRef);
    }
    m.matchScoreL = anchorL - getMismatchScore(headRef, m.caRead, NULL, 0, anchorL);
    m.mismatchScore = getMismatchScore(headRef, m.caRead, m.rawScores, 0, anchorL);
    m.matchScoreR= anchorL - getMismatchScore(&tailRef[examLength - anchorL], &(m.caRead[examLength]), NULL, 0, anchorL);
    m.mismatchScore += getMismatchScore(&tailRef[examLength - anchorL], &(m.caRead[examLength]), m.rawScores, 0, anchorL);

    if (m.mismatchScore <= (int)misT) {
        const char* ref1 = &(headRef[anchorL]);
        const char* ref2 = tailRef;
        int splicePos = (int)anchorL + findBestSplicedSite(ref1, ref2, m);
        if (m.mismatchScore <= (int)misT) {
            getSepCigarStr(m, splicePos, (int)anchorL);
            // Copy the middle part
            int examMidLength = (int)m.uiReadLength - 2 * anchorL;
            getSpliceAlignedRef(ref1, ref2, splicePos - anchorL, examMidLength, &(m.caRef[anchorL]));
            // Copy two anchor regions
            strncpy(m.caRef, headRef, anchorL);
            strncpy(&(m.caRef[examLength]), &tailRef[examMidLength], anchorL);
            m.caRef[m.uiReadLength] = '\0';

            m.MEMES = getMEMEScore4aMap(m, splicePos);
            if (m.MEMES >= this->opt.MEMES) {
                /*
                if (this->opt.bGetIntronSeq) {
                	string spliceSiteSeq = getSplicedSiteIntronSeq(ref1, ref2, splicePos, examLength);
                	strcpy(m.TAG, spliceSiteSeq.c_str());
                }*/
                bGoodDelMap = true;
            }
        }
    }
    return(bGoodDelMap);
}

inline bool checkStrEnd(char c)   // DEBUG
{
    if ( c == '\0') {
        cout << "End of a string" << endl;
        return(true);
    }
    return(false);
}

// This is a private funciton to get the intron seq in the border of junctions if available
// The intron sequence may be retrivable (if splice Pos is not exactly anchor regions)
// If not available, print *
string getSplicedSiteIntronSeq(string ref1, string ref2, int splicePos, int examLength)
{
    char IntronSeq[MAX_LINE];
    if (splicePos >= 2 && splicePos < examLength - 2) {
        IntronSeq[0] = ref1[splicePos + 1];
        IntronSeq[1] = ref1[splicePos + 2];
        IntronSeq[2] = ref2[splicePos - 2];
        IntronSeq[3] = ref2[splicePos - 1];
        IntronSeq[4] = '\0';
        // ref1[splicePos + 1], ref1[splicePos + 2], ref2[splicePos - 2], ref2[splicePos - 2]
        return(string(IntronSeq));
    } else {
        return(string("*"));
    }
}

int getSpliceAlignedRef\
(const char* ref1,  const char*  ref2, int splicePos, int refLength, char* caRefBuf)
{
    // splicePos is the length of the left arm - the left anchor region;
    int leftFragL = splicePos;
    int rightFragL = refLength - splicePos;
    strncpy(caRefBuf, ref1, leftFragL);
    strncpy(&(caRefBuf[leftFragL]), &(ref2[leftFragL]), rightFragL);
    caRefBuf[refLength] = '\0';
    return(0);
}

int getMEMEScore4aMap(CMappingResult& m, int splicePos)
{
    // splicePos is the length of the left arm - the left anchor region;
    int noMisMatchedL = 0;
    for (int i = 0; i < splicePos; i++) {
        if (m.caRead[i] != m.caRef[i]) {
            noMisMatchedL ++;
        }
    }
    int noMisMatchedR = 0;
    for (int i = splicePos; i < (int)m.uiReadLength; i++) {
        if (m.caRead[i] != m.caRef[i]) {
            if (m.caRead[i] != m.caRef[i]) {
                noMisMatchedR ++;
            }
        }
    }
    int leftArmL = splicePos;
    int rightArmL = (int)m.uiReadLength - splicePos;
    m.matchScoreL = getMatchScore(leftArmL, noMisMatchedL);
    m.matchScoreR = getMatchScore(rightArmL, noMisMatchedR);
    m.MEMES = min(m.matchScoreL, m.matchScoreR);
    return(m.MEMES);
}

inline bool checkCanonical(const char*  ref1,  const char*  ref2, int examLength, int splicePos)
{
    if (splicePos < 2 || splicePos > examLength - 2) {
        return(false);
    } else {
        char intron5_1 = ref1[splicePos];
        char intron5_2 = ref1[splicePos + 1];
        char intron3_1 = ref2[splicePos - 2];
        char intron3_2 = ref2[splicePos - 1];
        if (intron5_1 == 'G' && intron5_2 == 'T' && intron3_1 == 'A' && intron3_2 == 'G') {
            return (true);
        } else if (intron5_1 == 'C' && intron5_2 == 'T' && intron3_1 == 'A' && intron3_2 == 'C') {
            return (true);
        }
    }
    return(false);
}

// ref1 and ref2 are the two (middle) fragments that the read map to.
// m.uiPOS and m.uiMPOS are the left end indices of the two fragments.
// m.caRead should have the sequence; m.rawScores must should have the quality score
int CReadsDelMapping::findBestSplicedSite(const char*  ref1,  const char*  ref2, CMappingResult& m)
{
    int anchorL = (int)this->opt.anchorLength;
    int examLength = m.uiReadLength - 2 * anchorL;
    if (examLength < 0) {
        ERR
    }
    const char* readMiddleFrag = &(m.caRead[anchorL]);
    int splicePos = examLength;
    int minMismatchScore = getMismatchScore(ref1, readMiddleFrag, m.rawScores,  0, examLength);
    int misMatchScore = minMismatchScore;
    bool bCanoicalOnly = this->opt.bCanoicalOnly;
    bool bBestIsCanoical = false;
    for (int i = examLength - 1; i >= 0 ; i--) {
        if (minMismatchScore == 0)
            break;
        if ( checkStrEnd(ref1[i]) ||\
                checkStrEnd(ref2[i]) ||\
                checkStrEnd(readMiddleFrag[i]) ) {
            return(false);
        } // Assert

        if (ref1[i]  != readMiddleFrag[i]) {
            int baseScore = ( m.QScores[0] == '\0') ? 1 : (int)(m.rawScores[i]); //TODO check if the score is set
            misMatchScore -= baseScore;
        }
        if (ref2[i]  != readMiddleFrag[i]) {
            int baseScore = ( m.QScores[0] == '\0') ? 1 : (int)(m.rawScores[i]);
            misMatchScore += baseScore;
        }
        bool bCurrentRecordIsCanonical = checkCanonical(ref1, ref2, examLength, i);
        bool bBetterRecord = (misMatchScore < minMismatchScore);
        if (bCanoicalOnly) {
            if (bCurrentRecordIsCanonical) {
                if (!bBestIsCanoical) {
                    bBetterRecord = true;
                }
            } else {
                if (bBestIsCanoical) {
                    bBetterRecord = false;
                }
            }
        }
        if (misMatchScore == minMismatchScore) {
            bBetterRecord = !bBestIsCanoical && bCurrentRecordIsCanonical;
        }
        if (bBetterRecord) {
            splicePos = i;
            minMismatchScore = misMatchScore;
            bBestIsCanoical = checkCanonical(ref1, ref2, examLength, splicePos);
        }
        if (minMismatchScore == 0) {
            break;
        }
        // TODO: check more to see if there is a break of better MEMES or canoical
    }
    m.mismatchScore += minMismatchScore;
    return(splicePos);
}

// Calculate uiSkippedLength and CIGAR based on m.uiMPOS, m.uiPOS, m.strand and splicePos
int getSepCigarStr(CMappingResult& m, int splicePos, int anchorL)
{
    int armLength1 = splicePos;
    int armLength2 = m.uiReadLength - armLength1;
    //m.uiSkippedLength = m.uiMPOS - m.uiPOS - anchorL;
    if (m.strand == '+') {
        unsigned int leftEndOfSplice = m.uiPOS + armLength1;
        unsigned int rightEndOfSplice = m.uiMPOS + anchorL - armLength2;
        m.uiSkippedLength = rightEndOfSplice - leftEndOfSplice;
        sprintf(m.CIGAR, "%dM%uN%dM", armLength1, m.uiSkippedLength, armLength2);
    } else {
        unsigned int leftEndOfSplice = m.uiPOS + armLength2;
        unsigned int rightEndOfSplice = m.uiMPOS + anchorL - armLength1;
        m.uiSkippedLength = rightEndOfSplice - leftEndOfSplice;
        sprintf(m.CIGAR, "%dM%uN%dM", armLength2, m.uiSkippedLength, armLength1);
    }
    m.uiMPOS = 0;
    return(0);
}

inline bool validRefName(const char* RefName, int MAX_LENGTH)
{
    int i;
    for (i = 0; i < MAX_LENGTH; i++) {
        if (RefName[i] == '\0') {
            return(true);
        }
    }
    return(false);
} // For assert

// This private function extract the mapped info from the mapped info of two end
inline void copyMappedInfo(CMappingResult& m1, CMappingResult& m2, CMappingResult& m)
{
    m.uiRefLength = m1.uiRefLength;
    m.uiGlobalMappedPos = m1.uiGlobalMappedPos;
    m.uiPOS = m1.uiPOS;
    m.uiMPOS = m2.uiPOS;
    strcpy(m.RNAME, m1.RNAME);
}

int CReadsDelMapping::dealForwardDelAlignment \
(const CGenome_Index_TableQ& table, CMappingResult& m, CDelMappingsQ& delMappingQ)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];

    CMappingResult m1(q1, opt.anchorLength);
    CMappingResult m2(q2, opt.anchorLength);
    m.strand = '+';
    unsigned int hitNo1 = min(q1.ForwardAlignmentLoad, (unsigned int)q1.iMaxCapacity);
    for (unsigned int i = 0; i < hitNo1; i++) {
        if (i >= q1.iMaxCapacity) {
            ERR
            return(0);
        }
        getSingleMappingIndex(*table.pgenomeNT, q1, i, m1);
        unsigned int hitNo2 = min(q2.ForwardAlignmentLoad, (unsigned int)q2.iMaxCapacity);
        for (unsigned int j = 0; j < hitNo2; j++) {
            if (q1.aiHitIndex[i] < q2.aiHitIndex[j]) {
                getSingleMappingIndex(*table.pgenomeNT, q2, j, m2);
                if (strcmp(m1.RNAME, m2.RNAME) == 0  && m1.strand == m2.strand) {
                    copyMappedInfo(m1, m2, m);
                    if (this->examDelAlignment(*(table.pgenomeNTInBits), m, this->opt.subDiffThreshold)) {
                        // m.ISIZE is just a preliminary filtiering, due to we use paired-end setting for both end.
                        // The true separation should be in m.uiSkippedLength.
                        if (opt.disLB <= m.ISIZE && m.ISIZE <= opt.disUB) {
                            // update mapping result
                            delMappingQ.update(m, opt.bUniqueMapOnly, opt.bGetAllDelAlignments);
                            if (opt.bUniqueMapOnly &&\
                                    delMappingQ.bAmbiguousDelMapping &&\
                                    (delMappingQ.minDiff == 0)) {
                                return(2); // short cut
                            }
                        }
                    }
                } else {
                    break; // leave the for loop to another reference
                }
            }
        }
    }
    int noDelMapped = (int)delMappingQ.delAlignments.size();
    return(noDelMapped);
}

int CReadsDelMapping::dealReverseDelAlignment \
(const CGenome_Index_TableQ& table, CMappingResult& m, CDelMappingsQ& delMappingQ)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];

    CMappingResult m1(q1, opt.anchorLength);
    CMappingResult m2(q2, opt.anchorLength);
    m.strand = '-';
    unsigned int hitNo1 = min(q1.load, (unsigned int)q1.iMaxCapacity);
    for (unsigned int i = q1.ForwardAlignmentLoad; i < hitNo1; i++) {
        getSingleMappingIndex(*table.pgenomeNT, q1, i, m1);
        unsigned int hitNo2 = min(q2.load, (unsigned int)q2.iMaxCapacity);
        for (unsigned int j = q2.ForwardAlignmentLoad; j < hitNo2; j++) {
            getSingleMappingIndex(*table.pgenomeNT, q2, j, m2);
            if (strcmp(m1.RNAME, m2.RNAME) == 0 && m1.strand == m2.strand) {
                copyMappedInfo(m2, m1, m);
                if (this->examDelAlignment(*(table.pgenomeNTInBits), m, this->opt.subDiffThreshold)) {
                    // m.ISIZE is just a preliminary filtiering, due to we use paired-end setting for both end.
                    // The true separation should be in m.uiSkippedLength.
                    if (opt.disLB <= m.ISIZE && m.ISIZE <= opt.disUB) {
                        // update mapping result
                        delMappingQ.update(m, opt.bUniqueMapOnly, opt.bGetAllDelAlignments);
                        if (opt.bUniqueMapOnly &&\
                                delMappingQ.bAmbiguousDelMapping &&\
                                (delMappingQ.minDiff == 0)) {
                            return(1); // short cut
                        }
                    }
                }
            }
        }
    }
    int noDelMapped = (int)delMappingQ.delAlignments.size();
    return(noDelMapped);
}

int CReadsDelMapping::examDelAlignments4aRead\
(const CGenome_Index_TableQ& table, CMappingResult& m)
{
    bool bSepMore = false, bSepLess = false;
    sortAnchorByHitPos(this->alignmentsQ[0]);
    sortAnchorByHitPos(this->alignmentsQ[1]);
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];

    CDelMappingsQ delMappingQ;

    if (q1.ForwardAlignmentLoad && q2.ForwardAlignmentLoad) {
        // DEBUG
        // cout << q1.ForwardAlignmentLoad << "*" << q2.ForwardAlignmentLoad << endl;
        this->dealForwardDelAlignment(table, m, delMappingQ);
    }
    if ( (q1.ForwardAlignmentLoad < q1.load) &&
            (q2.ForwardAlignmentLoad < q2.load) ) {
        // cout << q1.load << "*" << q2.load << endl;
        this->dealReverseDelAlignment(table, m, delMappingQ);
    }

    int mapPosNo = (int)(delMappingQ.delAlignments.size());
    if (opt.bUniqueMapOnly && (mapPosNo == 1) && !(delMappingQ.bAmbiguousDelMapping)) {
        vector<CMappingResult>::iterator it = delMappingQ.delAlignments.begin();
        printDelMappingInSam(this->AlignResult, *it, opt.bPrintRef);
    } else {
        vector<CMappingResult>::iterator it = delMappingQ.delAlignments.begin();
        for (; it != delMappingQ.delAlignments.end(); it++) {
            printDelMappingInSam(this->AlignResult, *it, opt.bPrintRef);
        }
    }
    if (mapPosNo == 0) {
        this->bookKeepDelMapping(bSepMore, bSepLess);
        if (this->opt.bPrintUnMappedReads) {
            this->dealMissedRead(m);
        }
    } else {
        if (delMappingQ.bAmbiguousDelMapping) {
            this->noOfAmbiguousDelMapping++;
        }
        this->noDelMapping++;
        this->iMapDiffCount[delMappingQ.minDiff]++;
    }
    return(0);
}
/*
int CReadsDelMapping::dealMissedRead(const char* readName, CLongReadInBits& read, const char* QScores)
{
	char caBuf[MAX_LINE];
	char qScores[MAX_LINE];
	read.decode(caBuf);
	if(QScores != NULL) {
		unsigned int i;
		for(i = 0; i < read.uiSeqLength; i++) {
			qScores[i] = QScores[i];
		} qScores[i] = '\n';
		sprintf(this->MissReads->caBufp, "@%s\n%s\n+\nqScores\n", readName, caBuf);
		this->MissReads->fflush();
	} else {
		sprintf(this->MissReads->caBufp, ">%s\n%s\n", readName, caBuf);
		this->MissReads->fflush();
	}
	return(0);
}
*/
inline void set_ISIZE(CMappingResult &m1, CMappingResult &m2, int readLength)
{
    if (m1.strand == m2.strand) {
        m1.ISIZE = abs((int)m2.uiPOS - (int)m1.uiPOS);
    } else {
        m1.ISIZE = abs((int)m2.uiPOS - (int)m1.uiPOS) + readLength;
    }
    m2.ISIZE = m1.ISIZE;
}

void CReadsDelMapping::printInDelAlignStats\
(ostream& out, string readSetName, unsigned int uiSubThreshold)
{
    out << '\n';
    out << readSetName.c_str() << ", #Del-Mapping, " << this->iReadCounter << ", ";
    out << "#Filtered Reads, " << this->iBadReadCounter << ", ";
    out << "#Deletion Supported Reads, " << this->noDelMapping << ", ";
    out << "#Multi-aligned Reads, " << this->noOfAmbiguousDelMapping << "\n";
    out << readSetName << ", ";
    out << "#Deletion shorter than expectation, " << this->noOfDelSepMore << ", ";
    out << "#Deletion longer than expectation, " << this->noOfDelSepLess << ", ";
    out << "#Deletion shorter and longer than expected, " << this->noOfDelSepMoreAndLess << endl;
    out << readSetName << ", ";
}

// Note after the hit is sorted by the index, the q.asdiff is not the right value.
// This is ok because in the checking step, the whole reads will be evaluate again
void sortAnchorByHitPos(CAlignmentsQ& q)
{
    std::sort(&(q.aiHitIndex[0]), &(q.aiHitIndex[q.ForwardAlignmentLoad]));
    std::sort(&(q.aiHitIndex[q.ForwardAlignmentLoad]), &(q.aiHitIndex[q.load]));
}

// This read
void getLongReadSeqNQscores(CMappingResult &m,  CLongReadInBits& read, const char* QScores)
{
    m.strand = '+';
    read.decode(m.caRead, READ_BUFFER_SIZE);
    if (QScores == NULL) {
        m.QScores[0] = '\0';
    } else { // If quality score are available, get quality score
        trQScores(m.uiReadLength, SolexaScoreEncodingShift, QScores, m.QScores);
        m.getReverseReadandQual(); // TODO check this if it is correct
    }
}
